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Ab initio molecular dynamics calculations of deuterium desorbing from Si(lOO) 
have been performed in order to monitor the energy redistribution among the 
hydrogen and silicon degrees of freedom during the desorption process. The calcu- 
lations show that part of the potential energy at the transition state to desorption 
is transferred to the silicon lattice. The deuterium molecules leave the surface 
vibrationally hot and rotationally cold, in agreement with experiments; the mean 
kinetic energy, however, is larger than found in experiments. 



1 Introduction 

Hydrogen adsorption on and desorption from Si surfaces are of great technolog- 
ical relevance for, e.g., the etching and passivation of Si surfaces or the growth 
of Si crystals (see, e.g., Ref. [1] and references therein). Besides, the dynamics 
of the hydrogen interaction with Si surfaces is also of fundamental interest due, 
among others, to the so-called barrier puzzle: Mftile the sticking coefficient of 
molecular hydrogen on Si surfaces is very smallBBa indicating a high barrier to 
adsorption, in desorption experiments an almost thermal energy distribution 
of the molecules was found tl indicating a low barrier to adsorption. In order 
to explain this puzzle it was suggested to take tha^strong surface rearrange- 
ment of Si upon hydrogen adsorption into account: EfQ The hydrogen molecules 
impinging on the Si substrate from the gas phase encounter a Si configura- 
tion which is unfavorable for dissociation, while desorbing hydrogen molecules 
leave the surface from a rearranged Si configuration with a low barrier in the 
hydrogen coordinates. The lattice rearrangement energy was assumed to be 
about 0.7 eVB a value in reasonabLa^agreement with the transition state energy 
obtained by cluster calculationsElEltl Density-functional calculations of thejjit 
tential energy surface (PES) of H2/Si(100) using the supercell approach EHlEil 
provided a detailed microscopic description of the interaction. The transition 
state energy was determined to be only ^ 0.3 - 0.4 eV with the substrate 
rearrangement energy being merely < 50% of this value. 

The calculation of PESs is an important prerequisite for understanding 
reaction dynamics. For a quantitative analysis, however, a calculation of the 
dynamics is indispensable. We have therefore performed ah initio molecular 
dynamics calculations to monitor the energy distribution of D2 molecules des- 
orbing from Si (100). We will show that part of the potential energy at the 
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barrier position is indeed transferred to the silicon lattice, the mean kinetic 
energy of desorbing molecules, however, is larger than found in experiments. 

2 Calculational details 

In our ab initio molecular dynamics calculations the forces necessary to 
integrate the classical equations of motion are determined by density- functional 
calculations. The exchange-correlation functional is treated in the generalized 
gradient approximation (GGA) Ej. In previous slab studies the total energies 
were calculated within the local density approximation (LDA) with a posteriori 
GGA corrections Ej. The main effects of using the GGA in the complete seLL 
consistent cycle are a small increase of the theoretical lattice constant of SiliJ 
and a slight rise in the barrier height from Ei, — 0.3 eVtil to Eb — 0.4 eV. To 
correctly represent the up and down buckling of the clean Si(lOO) surface we 
use a (2x2) surface unit cell. The Si slab consists of five atomic layers. The 
topmost three of them are free to move in the molecular dynamics simulations, 
while the remaining two layers are fixed at their bulk positions. The density- 
functional calculations are performed with two k-points in the irreducible part 
of the Brillouin zone and 40 Ry cutoff energy. The equations of motion are 
numerically integrated within a predictor-corrector scheme with a time step of 
1.2 fs. The calculations have been performed on typically 64 nodes of the Cray 
T3D of the Konrad-Zuse-Zentrum, Berlin. 

3 Results 

Since the barrier tQ_associative desorption of hydrogen from Si (100) is rather 
high {Ed = 2.5 eV)Pj there is no sense in performing molecular dynamics cal- 
culations starting with the deuterium atoms at the adsorption sites because 
of the extremely low number of desorption events. Therefore we started the 
desorption trajectories close to the transiticpi. state for dissociative adsorption 
which was determined in the earlier studyO In total we have computed 40 
trajectories of D2 desorbing from Si (100). Eight trajectories were determined 
with the Si lattice initially at rest, i.e. at a surface temperature of = K, 
in order to specifically monitor desorption trajectories starting at the transi- 
tion state. Figure |^ illustrates the strong surface rearrangement of Si (100) 
upon hydrogen adsorption/desorption. Snapshots of a calculated trajectory 
are shown in Fig. ^). The dark Si atoms correspond to the relaxation of the 
Si lattice after the desorption event. Approximately 0.1 eV of the potential 
energy at the transition state is transferred to vibrations of the Si lattice which 
is a rather large amount compared to hydrogen/metal systems. 
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Figure 1: a) Hydrogen covered Si (100) surface (monohydride). b) Snapshots of a trajectory 
of D2 desorbing from Si (100) starting at the trasnsition state with the Si atoms initial at 
rest. The dark Si atoms correspond to the Si positions after the desorption event, c) Clean 

anti-buckled Si (100) surface. 



The desorption experiments determining the mean kinetic energy of D2 
desorbing from Si (100) were done with a rather high surface temperature of 
Ts ~ 920 K. In order to simulate these experimental conditions, we have per- 
formed ah initio molecular dynamics calculations with initial conditions corre- 
sponding to the experimental surface temperature. The system was allowed to 
equilibrate for more than 500 fs, whereby the deuterium atoms were kept close 
to the transition state by an additional cage-like potential. The additional po- 
tential was then switched off and the energy distribution of the D2 molecules 
desorbing from the thermal surface were monitored. In total 32 "thermal" 
desorption trajectories were calculated. The mean total, kinetic, vibrational, 
and rotational energies of the desorbing molecules are listed in Table |l| (note 
that UbTs = 0.079 eV). 

The results show vibrational heating, i.e. (Eyn,) > fc^Tsp-and rotational 
cooling, i.e. {Erot) < ksTs, in agreement with the experimentEj The mean ki- 
netic energy, however, is much larger than the experimental value of {E^inY^^ = 
0.165 eV. The difference between the experimental and theoretical results cor- 
resppnds roughly to the barrier height E\j. A closer analysis of the trajecto- 
ries Ej reveals that still approximately 0.1 eV of the potential energy at the 
transition state is transferred to the Si lattice, however, due to the Si lat- 
tice vibrations the mean adsorption barrier is increased by roughly the same 
amount. Possible contributions to the discrepancy between theory and exper- 
iment could be: (i) insufhcient statistics, i.e., too few trajectories computed, 
(ii) quantum mechanical effects (e.g., tunneling and zero-point effects) not 
taken into account in the classical molecular dynamics, (iii) dissipation chan- 
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Table 1: Mean energy distribution averaged over 32 trajectories of D2 molecules desorbing 
from a Si(lOO) surface at a surface temperature of Ts = 920 K. 



{Etot) 


0.72 ± 0.17 eV 


(Evib) 


0.11 ±0.09 eV 


{Ekin) 


0.58 ± 0.13 eV 


{Erot) 


0.03 ± 0.05 eV 



nels not considered, e.g., electronic excitations, (iv) limitations of the GGA 
functional, (v) experimental uncertainties. Certainly there is a strong need for 
future theoretical and experimental studies of the hydrogen on silicon system. 
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